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ABSTRACT 


The  performance  of  a  direction  of  arrival  estimation  procedure  at  low  signal- 
to-noise  ratios  and  limited  data  samples  is  an  important  characteristic.  The  approach 
based  on  maximum  likelihood  (ML)  estimation  is  considered  to  be  among  the  best  for 
this  problem  as  long  as  the  underlying  signal  model  is  properly  chosen.  Unfortunately,  in 
most  cases,  there  is  no  closed-form  solution  so  fast  search  procedures  are  employed.  Given 
no  a  priori  knowledge,  selecting  the  initial  parameter  values  for  these  search  procedures 
can  be  a  difficult  problem,  especially  under  low  signal-to-noise  conditions.  In  this  paper, 
a  new  method  for  uniform  linear  sensor  arrays  which  overcomes  the  initial  value  problem 
is  introduced.  This  method  is  called  Balanced  Data  Focusing  (BDF).  Simulation  results 
are  included  comparing  the  performance  of  this  new  method  to  that  of  the  ML  approach 
using  the  Alternating  Projection  Maximization  search  procedure  and  another  popular 
estimation  approach,  the  root-MUSIC  method. 

RESUME 

La  performance  des  techniques  permettant  I’estimation  de  Tangle  d’arrive  d’un 
signal  en  presence  d’un  faible  rapport  signal  sur  bruit  et  d’un  nombre  limite  d’echantillons 
est  une  characteristique  importante.  L’approche  se  fondant  sur  I’estimation  a  maximum  de 
vraisemblance  (ML)  est  consideree  comme  I’une  des  meilleurs  repondant  a  cette  exigence 
en  autant  que  le  mod^e  soit  adequatement  choisi.  Malheureusement,  dans  la  plupart  des 
cas,  une  solution  abregee  n’existe  pas  de  sorte  que  des  procedures  rapides  de  recherche  sont 
employees.  Sans  information  a  priori,  choisir  les  valeurs  initiales  des  parametres  pour  ces 
procedures  de  recherche  peuvent  s’averer  un  probleme  difficile,  specialement  dans  le  cas  ou 
on  a  a  faire  a  des  conditions  de  rapports  signal  sur  bruit  particulierement  bas.  Le  present 
rapport  introduit  une  nouvelle  methode  qui  s ’applique  a  un  reseau  lineaire  d’antennes 
uniformement  reparties.  Cette  methode,  appelee  BDF,  surmonte  le  probleme  du  choix 
des  valeurs  initiales.  Les  resultats  des  simulations  comprises  dans  ce  rapport  comparent 
les  performances  de  cette  nouvelle  methode  a  celles  de  deux  autres  methodes  populaires: 
la  methode  ML  employant  la  procedure  de  recherche  fondee  sur  la  maximization  des 
projections  alternees  ainsi  qu’a  la  methode  root-MUSIC. 
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EXECUTIVE  SUMMARY 


The  accuracy  of  a  direction  of  arrival  estimation  procedure  in  the  presence  of 
noise  is  an  important  characteristic,  especially  at  low  signal  to  noise  ratios.  Although  nu¬ 
merous  approaches  have  been  proposed  [1]  -  [5],  those  based  on  maximum  likelihood  (ML) 
estimation  are  considered  to  achieve  the  best  performance  [6]  assuming  the  underlying 
signal  model  is  properly  chosen. 

The  idea  behind  the  ML  approach  is  relatively  simple.  Using  an  appropriate 
model  of  the  signal  environment,  synthetic  data  can  be  generated  and  compared  to  the 
observed  data.  The  maximum  likelihood  parameter  estimates  are  taken  as  the  values 
of  the  model  parameters  for  which  the  synthetic  data  best  fits  (in  the  statistical  sense) 
the  observed  data.  Despite  the  simplicity  of  this  idea,  the  actual  derivation  of  the  ML 
approach  for  a  particular  problem  can  often  be  very  complex  depending  on  the  nature  of 
the  model  and  the  measured  data  available. 

For  the  bearing  estimation  problem,  the  deterministic  maximum  likelihood  (DML) 
method  has  been  shown  to  achieve  good  performance.  The  main  difiiculty  with  this 
method  is  that  it  requires  an  M  dimensional  search  where  M  is  the  number  of  signal 
bearings  to  be  estimated.  A  number  of  iterative  algorithms  have  been  developed  to  speed 
up  the  search  process  compared  to  a  brute  force  search.  In  general,  these  algorithms  ei¬ 
ther  suffer  from  a  requirement  for  good  initial  bearing  estimates  (without  which  they  may 
converge  to  the  wrong  solution),  or  their  convergence  time  is  still  too  slow  to  be  practical. 

To  overcome  these  problems  a  new  algorithm  has  been  developed  called  Balanced 
Data  Focusing  (BDF).  Rather  than  adjust  a  model  to  fit  the  data,  this  method  adjusts 
the  data  until  it  represents  a  valid  form  of  model  data.  The  bearings  can  then  be  easily 
calculated  from  the  data  using  a  linear  interpolation  method  (similar  to  linear  prediction 
methods).  The  advantage  is  that  initial  bearing  estimates  are  not  required,  and  this  is 
accomplished  without  sacrificing  convergence  speed  compared  to  other  approaches.  The 
main  restriction  is  that  the  sensor  array  is  required  to  be  linear  with  uniform  spacing 
between  sensors. 

In  simulations  the  BDF  algorithm  was  compared  to  the  deterministic  ML  method 
and  the  root-MUSIC  algorithm  (an  alternate  DF  estimation  approach  which  has  received 
considerable  attention  in  the  open  literature).  It  Wcis  found  that  the  BDF  method  per¬ 
formed  as  well  as  the  deterministic  ML  approach  and  outperformed  root-MUSIC  in  terms 
of  both  accuracy  and  threshold  performance.  The  BDF  method  is,  however,  slower  than 


V 


root-MUSIC,  especially  as  the  number  of  signals  increases.  Further  research  into  better 
search  algorithms  (compared  to  the  gradient  descent  technique  used)  would  likely  yield 
useful  improvements  in  processing  speed. 
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1.0  INTRODUCTION 


The  accuracy  of  a  direction  of  arrival  estimation  procedure  in  the  presence  of 
noise  is  an  important  characteristic,  especially  at  low  signal  to  noise  ratios.  Although  nu¬ 
merous  approaches  have  been  proposed  [1]  -  [5],  those  based  on  maximum  likelihood  (ML) 
estimation  are  considered  to  achieve  the  best  performance  [6]  assuming  the  underlying 
signal  model  is  properly  chosen. 

The  idea  behind  the  ML  approach  is  relatively  simple.  Using  an  appropriate 
model  of  the  signal  environment,  synthetic  data  can  be  generated  and  compared  to  the 
observed  data.  The  maximum  likelihood  parameter  estimates  are  taken  as  the  values 
of  the  model  parameters  for  which  the  synthetic  data  best  fits  (in  the  statistical  sense) 
the  observed  data.  Despite  the  simplicity  of  the  idea,  the  actual  derivation  of  the  ML 
approach  for  a  particular  problem  can  often  be  very  complex  depending  on  the  nature  of 
the  model  and  the  measured  data  available. 

In  this  report  a  linear  array  of  N  uniformly  spaced  narrowband  sensors  is  con¬ 
sidered.  The  complex  baseband  outputs  from  this  array  at  time  t  can  be  modelled  using 

M 

Vnit)  =  E  c„(t)e-^“-"^  (1) 

m=l 

where  M  represents  the  number  of  signals  impinging  on  the  array,  ?/„(<)  is  the  model 
output  for  sensor  n  (0  <  n  <  —  1),  Cm{t)  represents  the  complex  baseband  amplitude 

of  the  signal,  ujm  is  the  corresponding  spatial  frequency,  and  d  is  the  spacing  between 
sensor  elements.  The  spatial  frequency  is  related  to  the  signal  direction  of  arrival  <f)m  by 

27r 

(2) 

where  A  is  the  signal  wavelength.  The  direction  of  arrival  is  defined  as  the  angle  between 
the  sensor  baseline  (in  the  direction  from  sensor  0  to  sensor  A^  —  1)  and  the  direction  of 
the  signal  (which  is  perpendicular  to  the  wavefront).  The  model  is  adjusted  by  changing 
the  values  Cm{t)  and 

Thd  corresponding  measured  data  at  time  t  for  sensor  n  is  represented  by  Xn{t). 
This  data  will  differ  from  the  model  data  due  to  the  effects  of  noise,  mutual  coupling, 
multipath  propagation,  etc..  In  this  report,  only  noise  is  considered;  in  particular,  white 
Gaussian  noise,  since  this  is  the  most  common  type  of  noise  in  practice.  The  problem 


then,  is  to  estimate  the  bearings  of  the  incoming  signals  given  the  noisy  data. 

Using  these  definitions  of  the  model  and  measured  data,  various  DF  estimation 
methods  based  on  maximum  likelihood  theory  can  be  developed.  In  Section  2,  one  such 
method,  deterministic  maximum  likelihood  (DML),  is  discussed  along  with  the  difficulties 
associated  with  it.  This  is  followed  in  Section  3  by  a  description  of  a  proposed  new 
algorithm.  Balanced  Data  Focusing  (BDF),  which  overcomes  some  of  these  difficulties. 
Detailed  descriptions  of  various  aspects  of  the  BDF  algorithm  are  discussed  in  Sections 
4  and  5.  In  Section  6,  the  results  of  computer  simulations  are  used  to  compare  the 
performance  of  the  BDF  method  to  the  DML  method,  as  well  as  to  the  MUSIC  algorithm 
(an  algorithm  which  has  received  considerable  attention  in  the  open  literature).  The 
conclusions  are  presented  in  Section  7. 

2.0  THE  DETERMINISTIC  MAXIMUM  LIKELIHOOD  METHOD 

Assuming  that  has  the  same  form  as  yn{i)  with  the  addition  of  white 

Gaussian  noise  which  is  uncorrelated  from  sensor  to  sensor,  then  for  K  data  samples  the 
best  fit  between  the  model  and  data  is  achieved  when  [6]  the  mean  squared  error  (MSE) 
given  by 

A'-l  iV-l 

=  E  E  Mk)  -  2/n(4)|"  (3) 

A;=:0  n=0 

is  minimized.  This  is  called  the  deterministic  maximum  likelihood  (DML)  approach  since 
information  about  the  statistics  of  the  model  parameters  is  not  incorporated  into  the 
fitting  procedure  i.e.  the  statistics  are  assumed  to  be  unknown. 

For  the  purposes  of  this  reports,  it  is  also  convenient  to  represent  the  above 
expression  in  vector  form  as 

MSE  =  ^\x(4)  -  y{h))^Wk)  -  y{tk))  (4) 

k=0 

where  x{tk)  and  y  (4)  are  iV  x  1  elements  vectors  whose  elements  are  given  by  the  measured 
data  Xo{tk),Xi{tk),...,XN-i{tk)  and  the  model  data  yo(4),yi(4),  ■■■■,yN-i{tk)  respectively, 
and  the  superscript  H  represents  the  conjugate  transpose  operation.  The  vectors 
and  y(tk)  represent  samples  of  the  measured  and  modelled  data  taken  at  time,  tk,  and 
are  referred  to  as  snapshots. 

The  main  difficulty  with  determining  the  model  parameters  which  minimize  the 


2 


MSE  is  that  no  closed  form  solution  exists.  There  are,  however,  closed  form  solutions  for 
the  complex  amplitudes  (represented  by  Cm{tk))  given  that  the  values  for  the  bearings 
^2)  •  •  •)  [6]  are  known.  For  example,  the  signal  model  data  can  be  represented  as 

y(4)  =  Sc(4)  (5) 

where  S  is  the  N  x  M  normalized  signal  matrix  given  by 


^-jwi2d 


1, 

^—j(jj22d 


c(tk)  is  the  M  X  1  vector  of  complex  amplitudes  where 

Cl{tk) 

C2{tk) 

c(4)  =  03(4) 


CMitk) 


The  MSE  can  be  rewritten  as. 


MSE  =  Yj  “  Sc(4))^(x(4)  -  Sc(4)) 


Minimizing  this  expression  with  respect  to  c(tfc)  and  solving  leads  to  the  result  [6] 


c(4)  =  (S^S)-'S'")x(4) 


Plugging  this  result  back  into  equation  (8)  yields 


MSE  =  Y  -  S(S'^S)-iS^)x(4) 


where  I  is  the  N  x  N  identity  matrix. 

The  most  obvious  advantage  of  equation  (10)  is  that  it  simplifies  the  problem 
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from  a  search  for  the  complex  amplitudes  and  bearings  (3M  real  parameters)  to  a  search 
for  bearings  only  (M  real  parameters).  If  the  initial  estimate  of  the  bearings  can  be  made 
with  reasonable  accuracy,  several  good  methods  exist  which  can  fine-tune  the  estimates 
until  the  MSE  is  minimized  [6]  -  [9].  The  main  danger  of  these  methods  is  that  if  the  initial 
values  are  improperly  chosen,  they  will  converge  on  solutions  which  are  not  optimum  (i.e. 
a  local  minimum  of  the  MSE  function  instead  of  the  global  minimum). 

The  selection  of  initial  values  is  typically  performed  using  an  estimator  which 
requires  no  initial  values  but  is  less  accurate  than  the  DML  approach.  At  low  signal-to- 
noise  ratios  this  can  cause  problems  when  the  initial  value  estimates  become  so  inaccurate 
that  the  search  method  used  to  determine  the  DML  estimates  is  no  longer  guaranteed  to 
find  the  global  minimum.  The  results  are  therefore  unreliable  under  these  conditions. 

To  overcome  the  initial  value  problem,  an  exhaustive  search  of  all  possible  so¬ 
lutions  can  be  performed,  or  more  efficient  approaches  such  as  the  Simulated  Annealing 
algorithm  [10]  can  be  used  to  overcome  the  initial  Value  problem.  Unfortunately  these 
methods  achieve  better  reliability  by  sacrificing  convergence  speed  and  are  generally  too 
slow  for  most  practical  applications. 

3.0  THE  BALANCED  DATA  FOCUSING  ALGORITHM 

To  overcome  the  initial  value  problem  without  sacrificing  convergence  speed  a 
new  approach,  the  BDF  algorithm,  is  proposed.  To  simplify  the  following  discussion,  the 
single  snapshot  case  is  discussed  here  —  extensions  to  multiple  snapshots  are  presented 
in  Section  5.  Dropping  the  argument  tk  for  convenience,  the  equation  for  MSE  simplifies 
to 

MSE  =  (x  -  y)^(x  -  y)  (11) 

In  the  BDF  algorithm,  the  initial  value  problem  is  overcome  by  adjusting  or 
“focusing”  the  data  instead  of  adjusting  the  signal  model.  The  purpose  is  to  focus  the 
data  until  it  represent  a  valid  form  of  the  model  data,  i.e.  defining  S  as  the  set  containing 
all  possible  values  y,  the  focusing  procedure  continues  until  x'  e  S  where  x'  is  the  focused 
data. 

A  simple  approach  to  performing  the  focusing  is  to  use  a  gradient  descent  tech¬ 
nique.  Assuming  there  is  an  error  function  e(x)  that  provides  an  estimate  of  the  least 
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square  distance  from  x  to  S,  then  the  focusing  procedure  can  be  defined  by 

Xi+i  =  Xi  -  figi  (12) 

where  x,  represents  the  data  after  the  iteration,  Xq  =  x  the  initial  value  ,  //  the  step 
size,  and  g,  is  the  gradient  defined  by, 

gi  =  Vxe(x,)  (13) 

The  procedure  terminates  when  |e(x,)|  is  less  than  a  predetermined  tolerance  value. 

Since  the  set  S  contains  an  infinite  number  of  solutions,  there  will  be  an  infinite 
number  of  possible  solutions  for  the  focusing  procedure.  The  ideal  choice  is  the  solution 
which  lies  closest  to  Xo  in  the  least  squares  sense.  Therefore  the  development  of  the  error 
function  e(x)  is  critical  to  the  success  of  the  BDF  algorithm  and  is  discussed  in  Section 
3.1.  The  gradient  is  derived  in  Section  3.2.  Once  the  algorithm  has  converged,  any  number 
of  DF  algorithms  [l]-[5]  can  be  used  to  extract  the  bearings  from  x,.  However,  since  the 
calculation  of  e(x)  involves  the  estimation  of  the  autoregressive  filter  coefficients  of  the 
underlying  signal  model  in  x  (see  Section  4),  the  bearing  estimates  can  also  be  determined 
by  rooting  the  filter  coefficients  once  the  algorithm  has  converged.  The  step  size  fi  also 
plays  a  role  in  the  accuracy  and  convergence  time  of  the  algorithm  and  its  selection  is 
discussed  in  Section  3.3. 
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aM  = 

Equation  (14)  can  be  rewritten  as 


A/+l^j(wi+W2+W3  +  ...+WjVf)(i 


y~!  bmyn-M+m  0 


where 


bM  = 


It  is  relatively  straightforward  to  show  that  the  filter  coefficients  are  conjugate  symmetric, 


i.e.  bm  —  b\f_^- 


The  filter 


bi 

b  =  62 


is  a  spatial  moving  average  filter  with  M  nulls  corresponding  to  the  directions  of  the  M 
signals  in  the  model  data.  Applying  this  filter  to  the  measured  data  x,  the  signals  will  be 
nulled  out  and  the  resultant  output  will  be  noise  only.  That  is, 

M 

V  b*^Xn-m  =  Vn  (20) 


where  is  the  noise  or  error  output.  Applying  this  to  all  the  data  (i.e.  letting  n  = 
M,M  +  l,...,iV-  1)  then 

B^x  =  n  (21) 
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where 


0 


^0 

B  =  .  (22) 

^0 

0  6m 

and  n  is  an  (TV  —  M)  x  1  vector  with  elements  timiVm+i, 

Since  the  error  output  of  the  filter  will  be  a  function  of  the  signal-to-noise  ratio, 
then  an  appropriate  error  function  for  the  focusing  procedure  discussed  in  Section  2.0  is 
given  by  the  sum  of  the  square  errors.  Using  equation  (21)  this  can  be  expressed  as 

e(x)  =  n^n  =  x^BB^x  (23) 

The  filter  coefficients  6o,  6i, ...,  6m  could  be  generated  using  equations  (16),  (17), 
and  (18),  however  this  provides  no  advantage  since  initial  value  estimates  are  still  required 
for  the  bearings  <^i,  <^2,  •••5  <{>M-  Instead,  the  filter  coefficients  are  estimated  from  the  data 
itself.  The  methods  by  which  this  is  done  are  discussed  in  some  detail  in  Section  4. 

3.2  The  Gradient 

The  gradient  was  previously  defined  in  equation  (13).  Using  the  expression 
developed  for  e(x)  (equation  (23))  and  performing  the  gradient  operation. 

Si  =  Vx£(x£) 

=  Vx(xfBB^x.) 

=  2BB^X£  (24) 

3.3  The  Step  Size  fi 

An  appropriate  step  size  can  be  determined  using  a  first  order  Taylor  series 
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expansion  of  £(x,)  to  get 

£:(xi+i)  £(x,)  +  ~  2 

Substituting  the  right  side  of  equation  (12)  for  x^+i  and  solving  for  p  yields 

£(xi)  -  e(xi+i) 


To  minimize  convergence  time  e(xi+i)  =  0  is  chosen  which  gives 

,  ^  ^  (27) 

gf  g. 

Ideally,  with  this  choice  the  focusing  procedure  converges  in  one  step.  However,  since 
equation  (26)  is  only  an  approximation  a  smaller  step  size  is  required.  Accordingly, 

Po£(x,)  (Ofi\ 


where  0  <  po  <  1-  Through  empirical  testing  it  has  been  found  by  the  author  that 
po  =  0.3  provides  a  good  compromise  between  convergence  speed  and  accuracy. 

Using  this  last  result  the  expression  for  updating  the  data  (equation  (12))  be¬ 


comes 


Xi+i  =  Xi  - 


Po£(Xt)gi 


4.0  GENERATING  THE  FILTER  COEFFICIENTS 


For  the  purposes  of  this  analysis  Equation  (23)  can  be  rewritten  as 


<7^  =  b^XfX^b 
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where  cP"  =  e(x)  and  X/  is  the  forward  data  matrix  defined  by 


Xq, 

Xi, 

X2, 

•  5  XN-M-l 

Xi, 

X3, 

•  5  xn-m 

II 

X2-, 

X4, 

•  )  Xi^^M+1 

XMi 

Xm+1, 

Xm+2i  •  ■ 

• ,  xn^i 

The  goal  here  is  to  estimate  the  M  +  1  filter  coefficients,  6o,  6i,  62? --m  from  the  data 
instead  of  being  required  to  explicitly  specify  This  can  be  done  by  minimizing 

with  respect  to  b  subject  to  some  constraints. 

An  obvious  constraint  on  b  is  that  it  must  be  conjugate  symmetric,  since  fil¬ 
ters  generated  from  the  model  signals  have  this  property.  An  important  feature  of  the 
conjugate  symmetry  property  is  that  it  only  occurs  for 

(i)  complex  sinusoidal  signals  of  the  form  (see  also  equation  (1)) 

(32) 


(ii)  dampled  complex  sinusoidal  signal  pairs  of  the  form 

(33) 

or  combinations  of  the  two-signal  models  where  ajn{tk)  and  bm{tk)  are  the  complex  am¬ 
plitudes  and  p  is  the  real-valued  damping  factor.  The  conjugate  symmetry  constraint 
therefore  ensures  that  the  bearing  estimates  derived  from  b  will  be  based  on  solutions 
which  correspond  to  case  (i)  or  (ii).  Other  researchers  [11]  have  noted  that  the  inclusion 
of  case  (ii)  does  not  significantly  degrade  estimator  performance.  Results  supporting  this 
assertion  are  given  in  Section  5. 

Based  on  the  conjugate  symmetry  property  of  b,  equation  (30)  can  also  be 
written  as, 

cr"  =  b  (34) 
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where  Xt  is  the  backward  data  matrix  defined  by 


1 

to 

^N—Si 

H 

_ 1 

^iV-35 

•  ?  ^M-\ 

X(,  = 

^N-41 

I  . 

•  5  ^M-2 

(36) 

The  equivalence  of  equations  (30)  and  (35)  was  derived  assuming  b  was  conjugate  sym¬ 
metric,  however  this  equivalence  also  holds  for  any  choice  of  b  as  long  as  there  exists  some 
complex  constant  «  such  that  |q:|  =  1  and  the  vector  ab  is  conjugate  symmetric.  Since  the 
constant  multiplier  a  has  no  effect  on  the  filter  nulls  (and  therefore  the  estimated  signal 
bearings),  then  if  b  is  not  conjugate  symmetric  but  a  value  for  a  exists,  it  is  assumed 
that  b  is  multiplied  by  a  to  make  it  conjugate  symmetric. 

A  second  constraint  that  is  required  is  a  nontriviality  constraint.  The  purpose  is 
to  avoid  the  solution  6o  =  =  •••  =  =  0  which  is  obviously  meaningless.  Accordingly, 

two  such  constraints  have  been  investigated  in  this  report,  namely 

b^b  =  1  (36) 

which  leads  to  an  eigenanalysis  ba.sed  approach,  and 

b^Ufc  =  1  (37) 

which  leads  to  a  linear  prediction/interpolation  approach.  For  the  latter  approach  k 
represents  the  row  in  X/  or  Xj  to  be  predicted/interpolated  (0  <  A:  <  M)  and  Ufc  is  an 
M  -|- 1  column  vector  of  zeros  except  for  the  element  which  is  Uk  =  1. 

The  choice  of  the  nontriviality  constraint  influences  the  performance  of  the  BDF 
algorithm.  The  purpose  of  investigating  two  different  choices  is  to  illustrate  the  differences. 
Neither  choice  was  selected  based  on  optimality  criteria,  rather  they  were  selected  because 
they  are  commonly  used  constraints  in  DF  estimation  theory.  The  next  two  sections 
describe  how  these  constraints  are  incorporated  with  the  conjugate  symmetry  constraint 
into  a  solution  for  the  minimization  of  <r^. 
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4.1  Eigenvector  Approach 


For  the  eigen  analysis  approach,  advantage  is  taken  of  the  equivalence  between 
equations  (30)  and  (34)  to  recast  the  minimization  as 

2a'^  =  b"X/6Xf,b 

where  X/j,  the  forward-backward  data  matrix  is  defined  by 

X/6  =  [X/,X6]  (39) 

Minimizing  this  expression  subject  to  the  constraint  given  in  (36)  and  temporarily  ignoring 
the  conjugate  symmetry  constraint  leads  to  the  result  (also  known  as  Pisarenko  Harmonic 
Decomposition  [12]) 

XfbX%h  =  2a^h,  (40) 

where  vector  b  is  an  orthonormal  eigenvector  of  the  matrix  X/^X^  and  2a^  is  the  cor¬ 
responding  eigenvalue.  Since  the  solution  is  not  unique,  b  is  chosen  as  the  eigenvector 
corresponding  to  the  smallest  eigenvalue  of  Rm-  Eigen  analysis  or  singular  value  decom¬ 
position  techniques  can  be  used  to  determine  the  eigenvectors  and  eigenvalues. 

Although  no  obvious  conjugate  symmetry  constraint  has  been  applied,  an  eigen¬ 
vector  of  the  desired  form  corresponding  to  the  eigenvalue  will  always  exist.  For 
example,  the  matrix 

Rft  =  Xf,X%  (41) 

has  the  structure 

^00)  ^10?  ^205  •  •  •  1  ^M,0 

no,  ni,  ^211  •••,  ^M-1,0 

no,  ni,  ^22)  •••)  ^M-2,0 

^fb  =  .  .  . 

^M-1,0,  »'M-2,1,  •••,  ?'io 

^Af,0,  fM-lfi,  ^M-2,0,  •  •  •  ,  no 

The  matrix  R/t  has  the  properties  that  it  is  both  Hermitian  symmetric  (r,j  =  rjj  and 
persymmetric  (r,j  =  Assuming  that  b  is  a  valid  solution  to  equation  (40) 

but  is  not  conjugate  symmetric  (nor  can  it  be  made  conjugate  symmetry  through  the 
multiplication  of  a  complex  constant)  then  it  is  apparent  from  the  form  of  R/(,  that  b 


(38) 
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must  also  be  a  valid  solution  where 


The  new  vector  b  will  also  be  an  orthonormal  eigenvector  and  b^b  =  0.  Using  these  two 
eigenvector  solutions  it  will  always  be  possible  to  construct  a  new  eigenvector  solution  b' 
which  is  conjugate  symmetric,  for  example, 

b'  =  (b  +  b).  (44) 

4.2  The  Balance  Linear  Interpolator 

The  solution  to  the  equation  (30)  using  the  linear  prediction/interpolation  non¬ 
triviality  constraint  (equation  (36))  is  given  by 

XyXf  b  =  (45) 

where  =  1,  is  an  (M  1)  x  1  vector  with  elements  uq  =  u\  =  ...  =  um  =  0  except 
Um  =  I,  and  is  the  prediction  error  given  by 

(T2  =  b^X/Xfb  (46) 

Since  the  multiplication  of  the  filter  coefficients  by  a  constant  phase  term  has  no 
effect  on  the  location  of  the  filter  roots  or  prediction  error,  an  alternate  form  of  equation 
(45)  which  more  easily  allows  the  conjugate  symmetry  constraint  to  be  imposed  is  given 

by 

XyXfb  =  (47) 

where  c  is  an  arbitrary  complex  term  such  that  |c|  =  1,  and  now  bm  =  c. 

The  conjugate  symmetry  constraint  can  be  imposed  by  noting  that  if  b  has 
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conjugate  symmetry,  it  also  satisfies  the  relationship 

XfcXf  b  =  cV^UM-m  (48) 

where  is  the  backward  data  matrix  defined  by 

^7V-3?  •  •  •  f 

^N—2^  3?  ^iV— 4’  *  •  •  ?  ^M—1 

—  ^N—3’>  ^N—4'>  ^7V— 5?  *  *  *  ?  ^M—2 

^0 

Therefore  if  b  is  conjugate  symmetric  it  must  simultaneously  satisfy  the  both  equations 
(47)  and  (48).  This  can  be  forced  by  adding  the  two  expressions  to  get 

XyXf  b  +  X/Xf  b  =  R/ib  =  2a\cUm  +  C*nM-m)  (50) 

where  R/6  was  defined  previously  (equation  (42)).  In  this  new  expression,  c  is  no  longer 
arbitrary  as  will  be  seen  shortly. 

Before  discussing  the  actual  solution  of  b,  cr^,  and  c,  it  is  useful  to  define  the 

matrix 

Q  =  RJ.  (51) 

where  the  elements  of  Q  are  given  by  qij  for  i,j  =  0, 1, M.  It  has  been  assumed  here  that 
R/t  is  invertible.  The  matrix  Q  has  the  properties  that  it  is  both  Hermitian  symmetric 
(%  =  ^*ii)  persymmetric  {qij  =  qM-i^M-j)  since  R/j  also  has  these  properties  [12]. 
These  properties  will  be  used  to  make  simplifications  in  the  following  analysis. 

From  this  definition  the  filter  coefficients  can  be  expressed  as 


(49) 


b  =  2<7^Q(cU™  +  C*UM-m) 
Using  the  fact  that  bm  =  c  then  u|^b  =  c  and 

2  _  _ _ 

2u^Q(|c|2u„  +  (c*)2uM-m) 


(52) 

(53) 
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Although  single  snapshot  processing  can  be  useful  in  some  situations,  in  many 
real  world  systems  multiple  snapshot  data  is  available.  Utilizing  such  data  can  greatly 
enhance  estimator  performance.  The  extension  of  the  BDF  method  to  multiple  snapshots 
is  relatively  straightforward. 

Reintroducing  the  argument  tk  to  represent  a  quantity  based  on  the  sensor  data 
measured  at  time  instance  then  for  the  multiple  snapshot  case  the  sensor  data  x(to)^ 
x(U),  ‘ .  *,  x(<A'_i)  will  be  available  for  processing.  For  convenience  the  sensor  data  matrix 
is  defined  as 

X  =  [x(fo),  x(ti), . . . ,  x(fK-i)]  (60) 

and  the  error  function  is  then  represented  by  £(X).  The  forward-backward  data  matrix 
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in  this  case  is  also  redefined  as 


X/6  =  [Xf,{to),Xfb{h),...,Xf,{tK-i)]  (61) 

where  Xfb{tk)  is  the  forward-backward  data  matrix  formed  from  a  single  data  snapshot 
x(tk)  using  equations  (31),  (35),  and  (39). 

The  determination  of  the  filter  tap  weights  proceeds  in  the  same  manner  de¬ 
scribed  in  Section  4  using  the  new  definition  of  X/j.  The  gradient  vector  given  in  equation 
(24)  is  easily  modified  to  become  the  gradient  matrix  equation 

G  =  2B^BX  (62) 


Finally,  the  expression  for  updating  the  data  (equation  (29))  becomes 


Xj+i  —  X,'  -f 


^o£(X,)G,- 

trace{Gf^G,} 


From  this  derivation  it  can  be  seen  that,  compared  to  the  single  snapshot  case, 
increasing  the  number  of  snapshots  to  K  increases  the  amount  of  processing  by  a  factor  of 
K.  For  a  large  number  of  snapshots  the  processing  requirements  may  become  prohibitive. 
These  requirements  can  be  eased  somewhat  by  noting  that  the  MSE  for  the  DML  approach 
given  in  equation  (10)  can  be  rewritten  as 

MSE  =  trace{XX^(I  -  S(S^S)-^S")}  (64) 

By  inspection,  it  is  clear  that  the  data  matrix  X  can  be  replaced  by  an  equivalent  matrix 
Z  as  long  as 

ZZ^  =  XX^  (65) 

For  the  BDF  algorithm,  this  means  that  the  columns  of  the  matrix  Z  can  be  processed 
in  place  of  x(to),  x(fi),  . . .,  x(t/<:_i).  If  the  choice  of  Z  is  restricted  to  square  matrices, 
then  the  maximum  effective  number  of  snapshots  that  need  to  be  processed  by  the  BDF 
algorithm  is  limited  to  N,  which  results  in  a  considerable  savings  in  processing  when 
K>N. 

One  approach  to  computing  the  matrix  Z  is  to  perform  an  eigen  decomposition 
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on  XX^  to  get 

XX^  =  XiVivf,  (66) 

j=0 

where  the  eigenvalues  are  ordered  so  that  Aq  >  Ai  >  A2  >  . . .  >  A^v-i  and  the  corre¬ 
sponding  eigenvectors  Vq,  Vi,  V2, . . . ,  Viv_i  form  an  orthonormal  basis  set.  The  columns 
of  Z  are  then  given  by 

for  e  =  0,1,2, 1  (67) 


Although  there  are  much  more  efficient  approaches  to  generating  Z  (e.g.  Cholesky 
or  QR  decomposition),  the  eigen  decomposition  approach  provides  a  simple  method  to 
introduce  further  improvements.  In  other  research  the  observation  has  been  made  that 
the  eigenvectors  Vq,  Vi, ...,  Vm-i  are  mainly  a  function  of  the  signal(s)  present  (i.e.  sig¬ 
nal  subspace)  and  the  eigenvectors  vm,  vm-i-i,  •••,  vjv-i  (i-e.  noise  subspace)  are  mainly 
a  function  of  noise,  e.g.  [4]  and  [5].  Therefore  the  noise  subspace  eigenvectors  can  be 
ignored  and  the  columns  of  Z  computed  according  to 

Zi  =  J\i-  A,,  Vi  for  i  =  0,1,2,  ...,M- 1  (68) 


where  A,,  compensates  for  the  cumulative  effects  of  noise  in  the  main  diagonal  of  the 
matrix  XX^  and  is  estimated  using. 


A  secondary  advantage  of  rejecting  the  noise  subspace  eigenvectors  is  that  the  number  of 
input  vectors  is  further  reduced  from  N  to  M. 


6.0  SIMULATION  RESULTS 


To  investigate  the  performance  of  the  BDF  approach,  computer  simulations  were 
run  based  on  a  linear  N  element  sensor  array  with  uniform  half-wavelength  spacing.  Noise 
was  additive  white  Gaussian  noise  with  variance  at  each  sensor  output.  The  noise  was 
also  considered  to  be  isotropic  in  nature  and  uncorrelated  from  sensor  to  sensor.  The 
signal-to-noise  ratio  (SNR)  is  defined  here  as  the  ratio  of  a  signal  with  unity  power  (i.e. 


16 


Cm(i)  =  1  in  equation  (1))  to  the  noise  power,  or 

SNR  = -10  log{ri^)  dB.  (70) 

In  the  simulations  used  to  generate  the  following  results  all  signals  were  constant  ampli¬ 
tude  and  unity  power  (unless  otherwise  stated).  For  uncorrelated  signals,  the  signal  phase 
for  each  signal  was  independently  and  randomly  chosen  between  —  tt  and  tt  radians  for 
each  snapshot.  For  fully  correlated  signals,  signal  phases  were  randomly  chosen  between 
— TT  and  TT  radians  for  each  snapshot,  but  the  relative  phases  between  signals  was  kept 
constant.  Each  trial  consisted  of  processing  K  snapshots  of  sensor  data  to  generate  a 
single  estimate  of  the  signal  bearings.  Successive  trials  consisted  of  using  the  same  signal 
data  but  increasing  the  contribution  of  the  noise  data  in  1  dB  steps  until  the  entire  range 
of  interest  was  covered.  Then  new  signal  and  noise  data  was  generated  and  the  process 
repeated  L  times. 

For  comparison  purposes  three  different  methods  were  used  to  estimate  the  signal 
bearings  for  each  trial.  These  three  methods  were  the  BDF  method,  the  DML  method, 
and  root-MUSIC  [15].  The  root-MUSIC  method  was  included  since  it  has  been  extensively 
reported  in  the  open  literature  and  become  a  standard  by  which  other  estimators  are  often 
compared. 

For  the  BDF  method  the  value  of  the  step  adjustment  parameter  fio  in  equa¬ 
tion  (63)  was  chosen  to  be  0.3.  It  was  empirically  found  that  values  greater  than  this 
degraded  accuracy,  while  smaller  values  increased  the  computation  time  without  signifi¬ 
cantly  increasing  accuracy.  In  the  event  that  an  iteration  resulted  in  a  greater  error  (i.e., 
£(X.+0  >  £:(X{)),  fi  was  decreased  by  30%  and  the  iteration  repeated.  Iterations  in  the 
algorithm  were  stopped  when  the  value  of  e(X.i)  <  10  ^^trace{X^X},  which  typically 
occurred  between  50  and  200  iterations.  This  value  was  more  than  sufficient  to  ensure 
accuracy  and  could  probably  have  been  relaxed  in  most  cases. 

For  the  DML  method  an  exhaustive  search  was  performed  by  varying  each  angle 
in  1  degree  increments  to  locate  minimums  in  the  MSE  function,  followed  by  a  fine  0.1 
degree  search  around  each  minimum  to  determine  the  global  minimum.  The  Alternating 
Projection  Maximization  algorithm  [7]  was  then  used  to  further  fine-tune  the  results. 

For  the  root-Music  estimator  the  estimated  (p-f  l)  x  (p-f  1)  covariance  matrix 
was  computed  using  equations  (31),  (35),  (39),  and  (41)  where  the  parameter  p  was  used 
in  place  of  M.  For  multiple  snapshot  processing  with  uncorrelated  signals,  the  value  of 
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p  was  chosen  to  be  A/’  -  1  (instead  of  M  as  for  the  BDF  method).  For  single  snapshot 
processing,  or  fully  correlated  signals,  the  value  of  p  was  chosen  to  optimize  accuracy  and 
threshold  performance  for  each  particular  signal  environment  (i.e.  spatial  smoothing  was 
used). 


Accuracy  is  defined  here  as 


Accuracy  = 


M  K  ^ 

£  ^{Ki  - 


m=l  /=1 


(71) 


where  L  represents  the  number  of  trials  and  is  the  estimate  of  the  bearing  for 
the  trial. 


6.1  Comparing  Estimates  from  Different  Methods 

Figure  1  compares  the  performance  of  the  BDF  method  with  the  DML  and  root- 
MUSIC  methods  over  100  trials  with  single  snapshot  processing  {K  =  1),  N  =  S  sensors, 
a  signal-to- noise  ratio  of  30  dB,  and  two  signals  with  bearings  of  40  and  120  degrees.  For 
the  root-MUSIC  method  the  value  p  =  5  was  used  since  this  not  only  optimizes  accuracy 
but  maximizes  the  number  of  signal  bearings  that  can  be  estimated  by  this  method  using 
a  single  snapshot.  In  this  example  the  BDF  method  produces  nearly  identical  results  to 
the  DML  method.  The  root-MUSIC  estimates  are  different  and  slightly  poorer  (in  this 
example  by  a  factor  of  1.4). 


6.2  Damped  Solutions 

In  deriving  an  error  function  for  the  BDF  method,  the  signal  model  also  incor¬ 
porated  undesired  damped  sinusoid  pairs  of  the  form  represented  by  equation  (33).  The 
undesired  solutions  introduce  a  mechanism  for  error  in  the  BDF  algorithm.  For  example, 
against  two  signals  with  different  bearings  the  DML  method  will  estimate  two  distinct 
bearings.  In  the  BDF  method,  if  equation  (33)  is  chosen  as  the  best  model  for  the  data, 
then  the  estimated  bearings  for  the  two  signals  will  be  identical.  This  is  obviously  not 
correct. 


In  practice  the  selection  of  the  undesired  forms  are  far  less  common  in  the  so¬ 
lution  than  the  desired  signals,  and  when  they  do  occur,  do  not  significantly  degrade  the 
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Accuracy  (deg)  Accuracy  (deg)  Accuracy  (deg) 
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Table  1:  Damped  sinusoid  solutions 


BDF  Signal  Pair 

DML  Estimate 

P 

4> 

<^1 

(j>2 

0.862 

67.61 

67.58 

67.66 

0.873 

68.14 

68.00 

68.15 

0.886 

67.00 

66.96 

67.08 

0.895 

66.92 

66.88 

67.00 

0.914 

67.52 

67.47 

67.56 

0.937 

67.65 

67.60 

67.68 

0.940 

67.67 

67.66 

67.73 

0.946 

66.85 

66.81 

66.90 

0.955 

67.47 

67.41 

67.49 

0.963 

67.61 

67.58 

67.65 

results.  To  illustrate  this  a  simulation  experiment  was  run  with  two  signals  at  bearings  of 
65  and  70  degrees  and  a  signal-to-noise  ratio  of  23  dB.  Using  single  snapshot  data  {K  —  1) 
estimates  of  the  signal  bearings  were  made  using  both  the  DML  and  BDF  methods. 

Out  of  500  trials,  damped  sinusoid  pair  estimates  occurred  3  times  for  the  BDF 
method.  The  estimated  bearing  <j)  and  corresponding  damping  factor  (f>  for  the  first  10  of 
these  undesired  solutions  is  shown  in  Table  1.  The  corresponding  DML  estimates,  (f)i  and 
<^2,  are  also  shown  in  the  table. 

It  is  readily  apparent  from  the  results  shown  in  Table  1  that  when  pairs  of 
undesired  signal  types  are  selected  by  the  BDF  method,  the  DML  solution  has  almost 
converged  to  a  single  bearing.  The  result  is,  that  compared  to  the  DML  method,  the 
bearing  accuracy  of  the  BDF  method  is  not  significantly  degraded. 


6.3  Comparison  of  the  EV  and  BLI  Approaches 

Figure  2  shows  the  performance  of  the  BDF  method  under  difficult  conditions 
(8  sensors,  single  snapshot  processing,  and  signals  at  closely  space  bearings  of  69  and  70 
degrees)  for  two  different  approaches  to  filter  estimation,  namely,  the  eigenvector  approach 
(EV)  and  the  balanced  linear  interpolator  (BLI).  For  this  and  the  following  simulations 
L  =  500  trials  (see  equation  (71)).  At  high  signal-to-noise  ratios  both  approaches  perform 
equally  well.  At  lower  signal-to-noise  levels,  a  point  is  reached  where  performance  degrades 
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at  a  significantly  higher  rate  than  would  be  predicted  from  performance  at  higher  signal- 
to- noise  ratios.  Defining  this  as  the  threshold  point,  then  in  this  example,  threshold 
occurs  at  51  dB  for  the  EV  approach,  and  47  dB  for  the  BLI  approach.  In  comparing 
the  BLI  approach  to  the  EV  approach,  the  BLI  clearly  outperforms  the  EV  approach. 
This  is  a  fortuitous  result  since  the  BLI  approach  is  computationally  less  intensive.  Given 
these  advantages,  the  following  simulations  were  performed  using  the  balanced  linear 
interpolator  (as  described  in  Section  4.2)  to  generate  the  filter  coefficients. 


6.4  Single  Snapshot  Estimator  Performance 

Figures  3  and  4  show  the  accuracy  of  the  BDF,  DML,  and  root-MUSIC  methods 
for  a  single-signal  environment  and  a  two-signal  environment,  respectively,  using  N  =  S 
sensors.  Not  surprisingly,  all  three  estimators  perform  almost  equally  well  (within  1  dB) 
in  the  single-signal  case.  The  performance  advantage  of  the  DML  or  BDF  estimators 
compared  to  root-MUSIC  becomes  apparent  when  multiple  signal  environments  are  con¬ 
sidered.  In  Figure  4,  for  example,  the  DML  and  BDF  estimators  have  better  accuracy 
than  root-MUSIC  for  signal-to-noise  ratios  above  threshold.  The  values  of  p  used  for 
root-MUSIC  in  Figures  3  and  3  were  p  =  7  and  p  =  5  respectively. 

Figures  5  and  6  show  two  examples  of  single  snapshot  processing  with  N  —  8 
sensors  where  two  signals  are  closely  spaced  (p  =  5  for  root-MUSIC).  In  both  cases  the 
threshold  points  for  the  BDF  and  DML  methods  occur  at  signal-to-noise  ratios  up  to  6 
dB  lower  than  root-MUSIC.  The  apparently  superior  performance  of  the  BDF  method 
compared  to  the  DML  method  for  signal-to-noise  ratios  below  the  threshold  point  for  each 
estimator  is  not  important,  since  the  bearing  estimates  in  this  region  are  unreliable  and 
often  meaningless  (random).  Therefore,  there  is  relatively  little  to  be  gained  by  comparing 
the  accuracies  of  different  estimators  for  signal-to-noise  ratios  below  the  threshold  point. 
The  superior  performance  of  root-MUSIC  above  threshold  in  Figure  6  is  surprising  since 
using  a  value  of  p  <  TV  —  1  effectively  reduces  the  array  aperture.  This  usually  results  in 
poorer  accuracy,  as  in  Figure  5,  however  in  Figure  6,  the  opposite  is  true. 

Figures  7  and  8  show  two  more  examples  of  single  snapshot  processing  with 
TV  =  4  and  TV  =  16  sensors  respectively  (with  p  =  2  and  p  =  10  respectively  for  root- 
MUSIC).  In  these  two  examples  the  threshold  points  for  all  three  estimators  is  identical. 
Above  threshold  the  BDF  and  DML  methods  have  superior  accuracy  to  root-MUSIC 
equivalent  to  a  2  dB  increase  in  the  SNR. 
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Signal  to  Noise  Ratio  (dB) 


Figure  3:  Estimator  bearing  accuracy  versus  SNR  (<^  =  40  degrees,  N  =  8,  K  =  1) 


Figure  8:  Estimator  bearing  accuracy  versus  SNR  =  29  and  31  degrees,  N 


Figure  9  shows  one  final  single  snapshot  processing  exanaple  for  8  sensors  and 
three  signals  (p  =  5  for  root-MUSIC).  Note  that  the  DML  results  are  not  shown  here  due 
to  the  computational  difficulties  of  performing  an  exhaustive  three  dimensional  search. 
Analysis  of  estimator  performance  on  a  trial  by  trial  basis  for  SNR  =  47  dB  and  SNR 
=  48  dB  has  shown  that  for  this  example  threshold  performance  of  the  DML  and  BDF 
methods  are  within  1  dB.  Above  threshold  their  accuracies  are  the  same. 


6.5  Multiple  Snapshot  Estimator  Performance 

The  performance  of  all  three  estimators  in  the  last  five  examples  is  quite  poor. 
Better  results  can  be  achieved  if  more  snapshots  are  available  for  each  trial.  Figure  10 
shows  the  performance  of  the  BDF  method  using  K  =  100  snapshots  per  trial  for  various 
types  of  noise  preprocessing  (see  Section  5.0)  and  assuming  the  signals  are  uncorrelated. 
These  include  no  preprocessing,  full  processing  using  equation  (66),  and  partial  processing 
using  equation  (66)  but  including  all  eigenvectors  (i.e.  i  =  0, 1, 2, . . . ,  A'  —  1  which  is 
equivalent  to  subtracting  the  diagonal  noise  matrix  A,,!  from  R^v).  From  the  results, 
performance  with  full  noise  processing  is  the  best,  so  that  in  the  following  examples  full 
noise  processing  for  the  BDF  method  is  performed. 


6.5.1  Uncorrelated  Signals 

Figure  11  shows  the  same  signal  environment  as  in  Figure  10  but  comparing 
the  BDF,  DML,  and  root-MUSIC  approaches.  Again  both  the  DML  and  BDF  methods 
outperform  root-MUSIC,  despite  the  fact  that  under  these  conditions  root-MUSIC  can 
utilize  the  full  array  (i.e.  no  spatial  smoothing,  or  p  =  7). 

Figure  12  shows  the  same  signal  environment  as  in  Figure  6  except  where  K  = 
100  snapshots  were  used  for  each  trial.  The  relative  performance  of  each  of  the  estimators 
is  the  same  as  in  Figure  12. 

Figure  13  replots  the  performance  of  the  BDF  method  from  Figures  6  and  12 
over  the  same  range  of  SNR.  This  figure  illustrates  the  improvements  that  result  with 
multiple  snapshot  processing  compared  to  single  snapshot  processing.  In  this  example, 
increasing  the  number  of  snapshots  K  from  1  to  100  is  equivalent  to  increasing  the  SNR 
by  40  dB. 
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Accuracy  (deg) 


Figure  12:  Estimator  bearing  accuracy  versus  SNR  (<^  =  69  and  70  degrees,  iV  =  8, 
K  =  100) 
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Figure  13:  BDF  method  bearing  accuracy  versus  SNR  for  K  =  1  and  K  =  100  {(j>  =  45 
and  55  degrees,  =  4) 
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Figure  15:  Estimator  bearing  accuracy  versus  SNR  {(j)  —  45,  47,  50,  60,  and  120  degrees, 
N  =  10,K  =  10) 


Figure  15  is  a  five-signal  (uncorrelated)  example  comparing  the  BDF  method  to 
root-MUSIC  {p  =  7)  for  the  case  where  the  ii"  =  10  snapshots  and  iV  =  10  sensors  were 
used.  Again  due  to  the  computational  difficulties  of  implementing  the  DML  method  for 
more  than  two  signals,  its  performance  is  not  shown.  Analysis  of  DML  performance  on 
a  trial  by  trial  basis  for  SNR’s  between  38  and  40  dB,  and  using  various  initial  starting 
bearings,  has  determined  that  threshold  for  the  DML  method  occurs  at  39  dB  or  higher. 
Above  threshold  the  DML  and  BDF  methods  were  found  to  have  the  same  accuracy. 


6.5.2  Correlated  Signals 

Figures  14  and  Figure  16  show  examples  of  increasing  the  number  of  snapshots  for 
the  same  signal  environment  used  in  Figure  7  and  assuming  the  signals  are  fully  correlated. 
In  Figure  14  the  threshold  for  the  BDF  method  occurs  at  a  lower  SNR  than  either  the 
DML  or  root-MUSIC  methods.  In  Figure  16  the  threshold  for  the  BDF  and  root-MUSIC 
methods  occur  at  a  slightly  lower  SNR  than  the  DML  method.  The  superior  threshold 
performance  of  the  BDF  method  compared  to  the  DML  method  can  be  explained  by  the 
fact  that  the  BDF  method  appears  to  favour  signal  geometries  with  closely  spaced  signals. 
This  assertion  is  supported  by  that  the  fact  that  in  Figure  4  the  threshold  performance  of 
the  BDF  method  was  slightly  worse  when  the  signals  were  widely  spaced  {<f>  =  4  and  60 
degrees).  At  higher  signal- to- noise  ratios  both  the  BDF  and  DML  methods  have  identical 
accuracies  and  outperform  root-MUSIC  due  to  the  loss  of  effective  aperture  as  a  result 
of  using  spatial  smoothing.  Spatial  smoothing  was  required  for  the  root-MUSIC  method 
due  to  the  fact  that  the  signals  were  fully  correlated. 

The  advantage  of  increasing  the  number  of  snapshots  K  is  also  illustrated  in 
Figure  17  which  summarizes  the  results  for  the  BDF  method  from  Figures  7,  14,  and  16. 
Although  the  improvement  is  smaller  than  for  uncorrelated  signals,  it  is  still  significant 
with  an  equivalent  SNR  improvement  of  approximately  20  dB  for  the  case  where  K  =  100 
compared  to  the  case  where  K  =  1. 


6.6  Unequal  Signal  Amplitudes 

Figure  18  shows  one  final  simulation  example  where  the  signal  amplitudes  were 
1.0  and  0.1  corresponding  to  signal  bearings  of  95  and  105  degrees.  The  signals  were 
assumed  to  be  uncorrelated  and  p  =  5  was  used  for  the  root-MUSIC  method.  Again 
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Figure  17:  BDF  method  bearing  accuracy  versus  SNR  for  K  —  1,  K  =  10  and  K  =  100 
{(j)  =  45  and  55  degrees,  N  —  i) 
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Figure  18:  Estimator  bearing  accuracy  versus  SNR  for  unequal  signal  amplitudes  {<f)  —  95 
and  105  degrees,  N  =  b,  K  —  100) 


the  BDF  method  performs  as  well  as  the  DML  method  and  better  than  the  root-MUSIC 
method  in  terms  of  threshold.  Above  threshold  all  three  estimators  have  identical  accu¬ 
racies. 

6.7  Processing  Speed 

Although  no  claim  is  made  that  the  algorithms  tested  were  optimally  coded,  or 
for  that  matter,  that  the  most  appropriate  search  procedures  were  chosen,  some  mention 
of  relative  processing  times  is  in  order.  Not  surprisingly  root-MUSIC  was  the  fastest 
algorithm  under  all  conditions.  The  BDF  algorithm  was  found  to  be  about  an  order  of 
magnitude  slower  than  root-MUSIC  due  to  the  slower  convergence  of  the  BDF  algorithm 
(typically  between  50  and  200  iterations).  For  the  single  snapshot  case  the  convergence 
time  of  the  BDF  method  was  found  to  be  relatively  insensitive  to  the  signal  parameters 
including  the  number  of  signals,  signal  amplitudes,  and  bearing  spacing.  For  the  multiple 
snapshot  case  convergence  time  was  slower  for  an  increasing  number  of  signals  due  to 
the  increase  in  number  of  computations  (oc  M)  and  a  general  increase  in  the  number  of 
iterations  required. 

As  stated  earlier,  the  DML  method  was  implemented  using  an  exhaustive  search 
to  approximately  locate  the  global  minimum  of  the  DML  error  function  (i.e.  determine 
the  initial  bearing  estimates),  followed  by  a  fine  search  to  fine  tune  the  results.  Ignor¬ 
ing  the  estimation  of  the  initial  bearings  for  the  moment,  fine  tuning  was  achieved  using 
the  APL-ML  algorithm.  This  algorithm  was  chosen  as  it  was  found  to  be  computation¬ 
ally  the  fastest  of  the  algorithms  surveyed  ([6],  [7],  and  [8])  for  the  case  of  two  closely 
spaced  signals.  Generally,  once  properly  initialized,  the  APM-ML  algorithm  converged 
in  a  slightly  greater  amount  of  time  than  the  BDF  method  for  the  two-signal  case,  al¬ 
though  occasionally  it  was  many  times  slower.  The  APM-ML  algorithm  was  significantly 
slower  for  cases  involving  three  or  more  signals  especially  when  they  were  closely  spaced 
in  bearing.  Returning  to  the  problem  of  estimating  the  initial  bearings,  the  additional 
processing  required  adds  to  the  total  DML  method  processing  time  compared  to  simply 
fine  tuning  the  results.  In  the  case  where  an  exhaustive  search  is  used  (as  reported  here), 
the  additional  time  can  be  an  order  of  magnitude  or  more  greater  than  the  time  required 
for  fine  tuning.  An  alternate  procedure  for  determining  the  initial  bearings  which  is  sim¬ 
ple  and  fast  was  proposed  in  [7]  .  Although  this  procedure  was  reported  to  work  well  in 
the  multiple  snapshot  case,  for  the  single  snapshot  processing  examples  described  in  this 
paper,  this  procedure  was  not  as  effective. 
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7.0  CONCLUSIONS 


In  this  paper  a  bearing  estimation  method  for  uniform  linear  arrays  called  Bal¬ 
anced  Data  Focusing  was  described  which  was  shown  to  achieve  maximum  likelihood 
performance.  The  method  uses  an  error  function  based  on  a  linear  interpolation  filter 
which  estimates  the  minimum  mean  squared  error  (MSE)  of  the  deterministic  maximum 
likelihood  (DML)  method.  Using  this  error  function,  data  input  from  the  sensors  is  ad¬ 
justed,  or  focused,  until  the  minimum  MSE  is  0.  The  focused  data  can  then  be  treated 
as  noiseless  data  and  any  number  of  techniques  used  to  extract  the  signal  bearing  infor¬ 
mation. 

In  simulations  it  was  found  that  the  BDF  method  performed  as  well  as  the 
DML  method  and  outperformed  root-MUSIC  in  terms  of  both  accuracy  and  threshold 
performance,  even  under  extreme  conditions  where  signals  were  very  closely  spaced  in 
azimuth.  An  important  advantage  of  the  BDF  method  compared  to  root-MUSIC,  is  that 
the  BDF  method  can  deal  with  both  uncorrelated  and  fully  correlated  signals  (or  the 
equivalent  single  snapshot  problem).  An  important  advantage  compared  to  the  DML 
method  is  that  no  initial  values  are  required.  A  poor  choice  of  inital  values  for  DML  can 
lead  to  meaningless  results. 

A  disadvantage  of  the  BDF  method  is  that  it  is  slower  than  root-MUSIC,  espe¬ 
cially  as  the  number  of  signals  increases.  Further  research  into  better  search  algorithms 
(compared  to  the  gradient  descent  technique  used)  would  likely  yield  useful  improvements 
in  processing  speed. 
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